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Abstract 

We study the ionic distribution near a charged surface. A new method for performing Monte 
Carlo simulations in this geometry is discussed. A theory is then presented that allows us to accu- 
rately reproduce the density profiles obtained in the simulations. In the weak-coupling regime, a 
theory accounts for the ion-image interactions, leading to a modified Poisson-Boltzmann equation. 
When the correlations between the ions are significant, a strong-coupling theory is used to calcu- 
late the density profiles near the surface and a Poisson-Boltzmann equation with a renormalized 
boundary condition to account for the counterion distribution in the far-field. 
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I. INTRODUCTION 



Study of charged surfaces in electrolyte solutions is of fundamental importance, since 
these can model lamellar liquid crystals, clays, biological membranes, electrodes, etc. In- 
teresting phenomena such as like-charge attraction between similarly charged surfaces has 



been observed in the presence of multivalent counterions 
theoretical j4-15|, simulational 



3|. There has been a great 



171], and experimental jl8| effort to clarifying the 



behavior of double layers near charged surfaces. In many approaches the theories assume 
that the entire system is composed of the same dielectric material. This, however, is not 
very realistic since clays, colloidal particles, and hydrocarbon membranes, have dielectric 
constant significantly smaller than that of the surrounding aqueous medium^ The dielectric 
discontinuity across the interface results in polarization effects 



13 



17 



19 



-|2l| which 



can significantly affect the ionic distribution near the surface. In the present chapter, we 
present a simple theoretical approach which allows us to accurately predict the counterion 
distribution near a charged wall which separates two environments with different dielectric 
constants. We consider separately the weak and the strong coupling regimes. Monte Carlo 
simulations are also performed in order to test our theoretical predictions. 

II. MONTE CARLO SIMULATIONS 



The simulations of long-range interacting systems are much more difficult than of systems 
with short-range forces. The difficulty is that one can not arbitrarily cut off the long- 
range Coulomb potential by using periodic boundary conditions, as is the case of the usual 

periodic images 



221 . For systems 



Lennard- Jones fluids. Instead one needs to consider an infinite number o:: 
of the system and then sum over these using Ewald summation methods 
with a planar geometry, such as an infinite charged wall in contact with an electrolyte, 
there is an additional complication which comes from the broken translational symmetry. 
In this section, we describe an approach that allows us to simulate such systems taking into 
account the dielectric discontinuity at the interface. The Monte Carlo (MC) simulations 
are performed in the NVT ensemble. The system is located in the right-hand half of a 
rectangular simulation box of dimensions X Ly X L^, centered at the origin of coordinate 
system. A charged wall of surface charge density —a is located ai z = 0. Nc = int [aL^y/qa] 



2 



neutralizing counterions of charge aq and effective radius Tc, are confined to the region 
< z < where q is the proton charge and a is the ionic valence. The dielectric 

constants on the two sides of the wall are different, given by ec and e^, for 2; < and 
2; > 0, respectively. Note that the dielectric discontinuity results in the appearance of the 
image charges in the region —Lz/2 < z < which will be discussed later. The Ewald 



summation 



22| is used in order to calculate the electrostatic potentials between the ions in 



the periodic replicas of the simulation box. To account for the slab geometry we use the 



correction proposed by Yeh and Berkowitz 23|. The complete derivation of the electrostatic 
energy is presented in the appendix. 

III. THEORY: WEAK REGIME 

We first present a theory that accounts for the results of the MC simulations in the 
weak coupling limit, when the characteristic Coulomb interaction between the counterions 
is smaller than the thermal energy, F = a^q^ / e^dksT ^ 1, where d is the characteristic 
distance between the condensed counterions. Using aq/r^d"^ = a, the plasma parameter 



becomes T = a/ a^g'^vra/ ewksT. The Bjerrum length is defined as = (3q^ /e.^, and is 7.2 A, 
for water at room temperature. 

Before studying the ionic distribution near a charged wall, we first need to understand the 
role of electrostatic correlations and the induced charges when cr = 0. To this end we consider 
a symmetric a:a electrolyte at concentration Cg confined to infinite half-space. Fig. [TJ The 
work necessary to bring an ion from the bulk to a distance Zq from the (uncharged) surface 
which separates the two regions with the different dielectric constants, and ec, can be 
calculated in terms of the electrostatic Green's function 24 1. 



To account for the interionic correlations and induced surface charge, we use the linearized 
Poisson-Boltzmann (Debye-Hiickel) equation. For symmetry reasons it is convenient to work 
in cylindrical coordinate system. Suppose that an ion of charge q is located at Zq, see Fig. [H 
The electrostatic potential inside the regions 1 and 2 satisfies 

V^0(s, z) - z) = -^6is)6{z - Zq) , (1) 

while in the regions 3 and 4 it satisfies the Laplace equation, 

VV(s,^) = 0, (2) 



FIG. 1: Representation of an electrolyte in the region z > r^- 



where k = y/Sira'^XBCs is the inverse Debye length. 

Writing the potential as a Fourier transform, 4>{s,z) = (l/47r^) f^^dk e'*^'* (l>{k,z), we 
obtain the following equation for the regions 1 and 2, 

f til 



and for the regions 3 and 4, 

^) , (4) 



9z 

where we have used a Fourier representation of the delta function 

r" + CXD 



oo 



Since the electrostatic potential must remain finite in the limits z — )■ oo and z — oo, we 
obtain the following solutions for each region: 



pz 

(6) 



^4{k, z) = A^^e 



kz 



where p = ^/k^+K^. 

To calculate the integration constants, we use the conditions of continuity of the elec- 
trostatic potential, 4>3{k,z) = 4>4{k,z) at 2; = 0, 02(^,2;) = 03(fc,2;) at 2; = fc and 



')2{k, z) = z) at z = Zq, and of the normal components of the displacement field, 



d<j)4{k,z) 
d4>2{k,z) _ 



d<j>z{k,z) 
' dz 
d4>2{k,z) 

dz 
d4>i{k,z) 



, at 2; = , 
= , at 2; = Tc , 
= Anaq , a.t z = Zq 



(7) 



dz ^ dz 

The last equation has been obtained by integrating Eq. |3] across the singularity at Zg. 
The Fourier transform of the electrostatic potential in the region 2 is found to be 



where 



i{k,z) 



2naq 



-p{Zq~z) 



+ e 



-p(z+Zq-2rc) 



f2{k) 



(8) 



fi{k) = pcosh {kvc) — fcsinh {kvc) H — ^psinh (kvc) — 

—k cosh (/cfc) , 



(9) 



f2{k) = pcosh {kvc) + /Esinh {kvc) H — ^psinh (fcrc) + 



—k cosh [kvc 



(10) 



and the inverse Fourier transform is 



2tt 



dk kJo{ks)(f)2{k, z) 



(11) 



where Jo{ks) is the Bessel function of order 0. 

We are interested in calculating the potential felt by an ion, located at distance Zq from the 
interface. Subtracting the self-potential q/ew{zq — z), after performing the explicit integration 
of the first term in Eq. [HI we find 



'^poiyZq) 



-^ + ^/ dke 



~w Jo 



-2v{z,-r,) k h{k) 

V f2{k) 



(12) 



Performing the Giintelberg charging process 25|], we obtain the work necessary to bring an 



ion from the bulk to a distance Zq from the interface 



W,{Zq) 



2et„ Jo P f2{k) 



A very accurate approximation to the above expression is 



WapiZq) 



(13) 



(14) 



This approximate form is much more convenient for numerical implementation |26|, |27|, 
since it requires calculating only one integral to determine Wi{rc) at the beginning of the 
calculation. 

We now return to the problem of interest. The system now is an infinite dielectric wall 
of charge density —a, located at z = 0, and the neutralizing counterions of charge aq and 
radius Vc, confined to < ^ < Lz/2. The dielectric constants are ec and e^, for 2 < and 
z > 0, respectively. For F < 1 (weak coupling limit) the electrostatic potential and the ionic 
density profile can be determined from the solutions of the modified PB equation 

_0 , . . 4:71 



z = -- 



-(j5{z) + aqp{z)] 



(15) 



where the counterion density is given by 



a e' 



-aql3cl>{z)-l3Wap{z) 



(16) 



aq J^"^^ dz e-°'?'9?^(^)-/3W/aj,{^) 

The ionic correlations and the surface polarization are taken into account through the po- 
tential Wap{z), with K = y^SirXBacr /qLz- In Fig. |2l we compare our results with the MC 
simulations, for various dielectric constants. As can be seen, the agreement between the 
theory and the simulations is excellent. 
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FIG. 2: The symbols are the simulation data, while the lines represent the solutions of the modified 
PB equation (Eq. [T5]) . The surface charge density is o" = 6.25 x 10^^ q/-^"^ and the monovalent 
counterion radius is rc = 2 A. 
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IV. STRONG COUPLING REGIME 



When r > 1, the mean field theory — such as the PB equation — is not able to ac- 
curately predict the ionic density distribtuion, because of the strong correlations between 
the counterions. In the limit F ^ 1, the counterions form a quasi-two dimensional strongly 
correlated liquid near the wall js], [28I, with an approximately hexagonal geometry Con- 
sider one counterion. The electric fields produced by the others counterions of the double 
layer approximately cancel each other. The counterion then interacts predominantly with 
the wall and with the ionic image charges, see Fig. [31 The potential produced by the charged 
plate which separates the two environments with different dielectric constants is given by 

, / N 47rcr ^ ^ 

M^) = -7-— 7T^ • (17) 

As an approximation, we consider that the ion interacts only with the self-image and with 
the image charges of the 6 first neighbors in the hexagonal lattice, see Fig. [3l 




aq ) (aq) ( aq 
aq) (aq 



FIG. 3: Hexagon of images at the surface. In (A) the side view. In (B) the self image and the 
nearest neighbors. In order to illustrate we consider ec = 0. 



30|. The 



This approximation was used previously in the study of colloidal double layers 
electrostatic energy of a counterion at distance z from the plate is then 

^,.) . ^ 2g , , (18) 

where 7 = (e^ — Cc) / (e^ + Cc) and h is the distance between the ions of the hexagonal lattice. 
h can be calculated by considering that = aA/aq ions are distributed on the surface of 



7 



area A. The unitary cell of a hexagonal lattice is a parallelogram of area /i^ a/3/2, which 
gives the result 

^ (19) 



h 



The ionic density profile near the surface is obtained from 



(20) 



where C = a/aqj^dz e"^^^^-* is the normalization constant. In Fig. H] we compare our 
theoretical results with the MC simulations. The agreement is very good in the region 
where the strong coupling approximation applies. 
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FIG. 4: The symbols are simulation data, while the lines represent the theory. The surface charge 
density is a = 3.74 x 10^'^ q/-^'^ and the pentavalent counterions radius is Tc = 2 A. The solid line 
in the inset shows the solution of the regular PB equation with the boundary condition given by 
Eq.EH 

In the far field, we expect that the counterions will be very dilute so that the electro- 
static potential will, once again, satisfy the PB equation. The boundary condition at the 
colloidal surface, however, must be modified to account for the strong counterion condensa- 
tion induced by the electrostatic correlations. The new boundary conditon can be derived 
by equating the electrochemical potential of the condensed counterions and of the counte- 



rions which remain in the bulk 



28 



31 



32| . This results in a new boundary condition for 



the standard PB equation which requires that the concentration of the counterions near the 
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surface be 

Ppb(O) = psce^''' , (21) 

where /3/ic = — 1.65r + 2.Q1T^^^ — 0.26 In F — 1.95 is the chemical potential of the strongly 
correlated counterions [29]. The density psc is obtained using the coarse-graining of the 
near- field density profile, Eq. [201 the region near the surface 3l| . 

where Xqc = l/27raABcr is the Gouy-Chapman length. In the inset of the Fig. IHwe present 
the solution of the usual PB equations with the renormalized boundary condition given by 
Eq. |2T1 Only the case with ec = is shown, since in the far field the ionic density distribution 
is highly insensitive to the value of ec- 



V. CONCLUSIONS 



We have presented a method for performing MC simulations in a cell geometry that 
includes a dielectric discontinuity at one of the boundaries. The results of the simulation 
have been used to study the counterion density profiles and to develop the weak and the 
strong coupling theories which account very accurately for the simulation data. In the weak 
coupling regime, the image charges repel the counterions from the wall. The contact density 
predicted by the present theory is substantially smaller than is found using the usual PB 
equation, and is in excellent agreement with the MC simulations. In the strong coupling 
limit, the contact density is found to be even lower, since in this case the counterions are 
repelled both by the self-image and by the images of the others counterions. Finally, we 
show how for F ^ 1 the counterion density distribution can be calculated in the far field 
using a renormalized boundary condition for the standard PB equation. 

In presenting the theory we have restricted ourselves to the systems containing only 
counterions and no colons. In the weak coupling limit, the approach developed here can be 
easily extended to systems which also contain 1:1 electrolyte. The situation, however, is much 
more difficult for multivalent electrolytes. For such systems, strong electrostatic interactions 
between the counterions and colons lead to formation of Bjerrum clusters. Thus, to be able 
to account for the distribution of multivalent ions near a charged surface one must first 
have an accurate description of the bulk of solution. This, already presents a formidable 

9 



challenge, see Ref. [3|. Nevertheless, one can make some progress by considering a chemical 
picture of electrolyte in which there is an equilibrium between the free ions and the clusters, 
such calculations, however, very rapidly become quite involved {32]. 
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Appendix A: Energy Calculation for Monte Carlo Simulations 

We consider a charge neutral system of ions of charges gj. The electrostatic potential 
at the position r, created by all ions (excluding ion i), their image charges (including the 
image of ion i), and the periodic replicas is 

00 N „ , s 



n j=l 



I e^\r-s + rep\ 



00 Af „ / 



n j=l 

where pj{s) = qj6{s — rj — r^p) and p'j{s) = 7gj5(s — r'j — r^p) are the charge densities of 
ions and their replicas; and of dielectric images and their replicas. The replication vector is 
defined as r^p = L^yUxX + L^yUyy + L^rizZ and r'j = rj — 2zjZ. The vectors n = {n^, Uy, Uz), 
where rix, riy and are integers, represent the infinite replicas of the main cell. The constant 
7 is defined as 7 = (e^ — ec)/(e«) + ^c) and the prime on the summation means that j 7^ i, 
when n = (0, 0, 0). The total electrostatic energy of the system is given by 

1 ^ 

U = -J2q^Mr^) . (A2) 

1=1 

The energy above is very difficult to calculate because of the slow convergence of the series 
in Eq. lAll To speed up the convergence, we use the Ewald method in which the ionic charge 
is partially screened by placing a Gaussian-distributed charge of opposite sign on top of each 



ion |22|. We then add and subtract opposite Gaussian charge at the position of each ion 
and its image, Pj{s) and p'j{s), respectively. The potential, Eq. lAlt then becomes 

Mr) = <Pf{r) + <P\r)-<pf^{r) , (A3) 
10 



where 



oo N 

EE' 

n j=l 



e,„ r -s + r 



ep I 



oo N 

EE 

n i=l 



PUS J 



eJr - s + r 



(A4) 



ep I 



e^„|r -s + r 



oo N „ iG / \ 



yy / , — .d^s (A5) 



and 



n j=l 



r) = / 7 , (A6) 

e,„r-s 



where pf(s) = (k^/v^) exp (-k^|5 - r^- - rgpp), Pj*^(5) = 

7Q'j('^e/V^) 6xp (— Kgis — r^- — TepP) and is a dumping parameter. We subtracted 
the self potential, Eq. IA6t from the Eq. IA3t in order to remove the prime over the 
summation in the long-range (L) part of the potential, Eq. IA5I The electrostatic potential 
produced by the Gaussian charges can be easily calculated using the Poisson equation, 
yielding 

^V)-EE^ /1if;V/T" + 



n j=l 

oo 



ep I 



2^ 2^ , \r-r' +r I ' ^^^^ 

where erf(x) is the error function. The short-range part of the potential (5*), Eq. IA4t can 
then be obtained in terms of the complementary error function, erfc(x) = 1 — erf(x), 

oo N 



lS^ \ V^V^' erfc(Ke|r -r- +r 

n j=l '"I ' J ^ ' ep 

oo N 



EE^^ 7 ;;;t" - m 

n j=l ^"'l 'j^'epl 

This potential decays very rapidly and can be truncated by setting the dumping parameter 
to Ke = where V = L'^yLz, corresponding to the minimum image convention. It is 
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then sufficient to consider in tlie sum only tlie term n = (0,0,0), witli tlie usual periodic 
boundary condition. 

The self-potential, Eq. IA61 reduces to 



erf(Ke|r -Til) 



(AlO) 



We next calculate the long-range part of the potential, Eq. IA7I This is most easily ob- 
tained using the Fourier representation, = (l/V) Jyd^r exp {—ik ■r)(f)^(r), since in 
the reciprocal space all the sums, once again, converge very rapidly. The Fourier transform 
p^{k) = (l/V) fy d?r exp {—ik ■ r)p^(r), of the Gaussian charge density, 



oo N 



P^ir) = V Vg^— ^exp(-K2|r-rj -r^pH + 

n j=l V /I 

oo N 3 



(All) 



IS 



n i=l 



1 IfcP 

p^(fc) = 77 exp (--^; 



^qjexp{-ik-Tj) + 
^7gjexp(-ifc-r^.; 



(A12) 



where = i2Tmx/ L^yj^T^ny/ Lxy,2T:nz/ L^). Using the Poisson equation, \k\ (j) (k) 
{4:TT / e^) (k) , we can evaluate the Fourier transform of the potential. 



\k) 



Air 



\k\ 



eu.V\k\' 



exp 



N 



qj exp {~ik ■ Tj 



N 



(A13) 



2^7gjexp(-ifc-rj.; 

i=i 

The corresponding real-space electrostatic potential is calculated using the inverse Fourier 
transform, 0^(r) = 0^(fc) exp {ik ■ r), 

^'^l'^) = exp (-x^) exp {ik ■ r) x 



tu,V\k\^ 



N N 

^ qj exp (-zfc ■ r^) + ^ exp (-ifc ■ r'-) 
Li=i i=i 



(A14) 
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The long-range contribution to the total electrostatic energy is given by Ul = 
(1/2) gj0-^(rj), where 4>^(r) is obtained from Eq. IA14I It is convenient to rewrite 
this in terms of functions: A{k) = Xlili 1i ^^^^ {k ■ ri), B{k) = — Xlili 1i ' ^i); C!{k) = 
cos (A; ■ and D{k) = —^^^i'yqisin{k-r'^. The electrostatic energy then be- 
comes, 



rr 27r 



2 

X 



k 

[A{ky + B{kY + A{k)Cik) + B{k)D{k)] . (A15) 

These functions are easily updated for each new configuration in a Monte Carlo sim- 
ulation. The electrostatic energy coming from the short-range part of the potential is 
Us = {^/'^)J2iLi1i'Pi(j'i): where 4>i(r) is given by the Eq. IA91 and the self-energy con- 
tribution is Useif = (1/2) Xlili Q'i^r''^ ('"«)• III the limit x — 0, the erf(x) function vanishes 
as (2/-\/7r)a; and the self-energy contribution reduces to, Ugeif = {'^e/^w\^)J2iLi1i- "^^^ 
total electrostatic interaction energy of the ions is given by the above expressions plus the 
correction for the slab geometry. Yeh and Berkowitz (23| found that the regular 3D Ewald 
summation method with an energy correction, can reproduce the same results as the 2D 
Ewald method, with a significant gain in performance. Taking into account the dielectric 
discontinuity and the induced image charges, we find the correction for the slab geometry 
to be 

" N 

-i=i 

N 

Y,l<lM-z',f , (A16) 

i=i 

where z'^ = —zj. Using the electroneutrality, this expression can be written as 

[/,„, = ^M,2(l-7) , (A17) 

where = Xlili li^i magnetization in the z direction. 

Now suppose that the system consists of Nc counterions of charge aq and a wall of uniform 
surface charge density —a, located at z = 0. We first derive the functions A, B, C and D 



TT ^ 
Ucor = — — TF ^ Q; 

^""^ i=l 
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appearing in the long-range part of the potential, Eq. IA15I For the surface charge we find 

/Lxy/'i I'Lxy/'i 
I a dx dy cos {k^x + kyu) = 

- -r-r^ sin {k^L^y/ 2) sin {kyL^^y/ 2) , 

/Lxy/2 pLxy/2 
I a dx dy sin (k^x + kyy) = , 

/Lxy/2 pLxy/2 
I 70" dx dy cos {k^x + kyy) = 
'Lxyl'^ J Lxy /'^ 

sin {kxLxy/2) sin {kyLxy/2) 



kxky 



and 



/Lxy/2 pLxy/^ 
/ 70" dx dy sin {k^x + kyy) = . 
~Lxy/'2 J Lxy/'2 

The corresponding functions for A^^ counterions and the charged wall are then: A{k) = 
agE£iCos(fc-r,) + Ap(fc), B{k) = sin (A; -r,), = 7«g E£ cos (fc ■ + 

Cp{k) and D{k) = —jaq sin {k ■ r-), and the total long-range part of the energy, Ul, is 
given by the Eq. RTSl 

The short-range contribution to the electrostatic potential created by the charged surface 
at distance Zi is 

(Pp{Zi) = --, ■ r X 



■Lxy/2 J~Lxy/2 + ?/2 + 

The limits of integration are defined in order to keep the minimum image convention. We 
calculate the potential on a grid in the z direction with spacing between the points 0.01 A. 
The calculation is performed once at the beginning of the simulation, and the potential is 
tabulated. The total short range electrostatic interaction energy is then given by Us = 
{aq/2) J^^^i (l)firi) + aqj^^^i M^i)^ where 0f (r) is 

^, erfcK|r-r,|) , ^ erfc(/.e|r - r;.|) 
ag> , i^ + iaq^ ■ . A19) 



The self energy can be written as Useij = ('^e/e^u^A")(^ctt^q'^ + cr^L'^y)- Since the charged 
surface is located at z = 0, it does not contribute to the correction potential, Eq. IA17[ so 
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that the magnetization remains — Oiq^-^^ Zi. The total energy used in the simulations 
is 

U^Us + Ul- Usei f + U^. (A20) 

We use 1 X 10^ MC steps to equihbrate the system. The configurations are saved each 100 
MC steps. The counterionic density profiles are obtained with 80 x 10^ saved uncorrelated 
states. 
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